% Evaluate total compliance by the 3-el WK method 

global time pressure flow T

colordef black;
aantal=length(pressure);

meanP=mean(pressure);
meanQ=mean(flow);

TVR=meanP/meanQ;
CO=meanQ*60/1000;
T=max(time);
HR=60/T;

% compute input impedance

Zin=fft(pressure)./fft(flow);
Zinmod=abs(Zin);
Zinphase=unwrap(angle(Zin));
freq=(1/T)*((1:aantal)-1);


% compute characteristic impedance in the frequency domain (average 3-15 Hz)

harm1st=1+ceil(3*T);
harmlst=1+floor(15*T);
Z0=mean(Zinmod(harm1st:harmlst));


% THE 3-element WK method

options=[];
options(14)=5000;
Rinit=TVR;
Z0init=Z0;
Cinit=1.5;

%est = fmins('fc_wk3',[Rinit Z0init Cinit]',options,[],time,pressure,flow,T);
est = fminsearch('fc_wk3',[Rinit Z0init Cinit]');

R_WK3=est(1);
Z0_WK3=est(2);
C_WK3=est(3);


% compute pressure respons and display ...

% step 1: transform Matlab fft flow into harmonics

fftQ=fft(flow);

fftQ(1)=fftQ(1)/aantal;

for j=2:aantal/2
    fftQ(j)=2*fftQ(j)/aantal;
end

% step 2: compute input impedance 3-el WK model and harmonics P_WK3 

for j=1:aantal/2
    freq(j)=(j-1)/T;
    Z_WK3(j)=Z0_WK3+R_WK3/(1+i*2*pi*freq(j)*R_WK3*C_WK3);
    fftPest(j)=Z_WK3(j)*fftQ(j);
end

% step 3: compute matlab harmonics

fftPest(1)=fftPest(1)*aantal;

% check for even or uneven number

if rem(aantal,2)==0  
	for j=1:aantal/2-1
	    fftPest(j+1)=fftPest(j+1)*aantal/2;
	    fftPest(aantal-j+1)=conj(fftPest(j+1));
	end;
	fftPest(aantal/2)=fftPest(j+1)*aantal/2;
end;
if rem(aantal,2)==1  
	for j=1:ceil(aantal/2)-1
	    fftPest(j+1)=fftPest(j+1)*aantal/2;
	    fftPest(aantal-j+1)=conj(fftPest(j+1));
	end;
end;


% step 4: compute inverse matlab transform

Pest=real(ifft(fftPest));

fig1=figure;
subplot(2,1,1);
semilogy(freq(1:11),abs(Zin(1:11)),'m+',freq(1:11),abs(Z_WK3(1:11)),'y');
hold on;
title('input impedance');
xlabel('frequency (Hz)');
ylabel('impedance (mmHg/(ml/s))');
legend('meas','model');

subplot(2,1,2);
plot(time,pressure,time,Pest);
hold on;
title('pressure');
xlabel('time (s)');
ylabel('pressure (mmHg)');
legend('meas','model');
fprintf(1,'\nestimated R: %g mmHg/(ml/s) Z0: %g mmHg/(ml/s) 3-elWK C: %g ml/mmHg\n',est(1),est(2),est(3))

